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MASSACHUSETTS 


ABSTRACT 


The  angle-of -arrival  estimation  problem  for  waves  incident  upon  a  sensor 
array  was  examined  through  a  Monte  Carlo  evaluation  of  the  performance  of  the 
Bayes-optimal  MAP  (maximum  aposteriori)  and  MMSE  (minimum  mean  square  error) 
estimators.  The  case  of  two  independent  wave  emitters  of  known  powers  as  well 
as  a  multiple  look,  Gaussian  signal  in  Gaussian  noise  statistical  model  were 
assumed.  The  Cramer-Rao  bound  on  the  estimator's  rms  error  was  computed  for 
comparison. 

The  evaluation  proceeded  with  the  computation  of  MAP  and  MMSE  angle  esti¬ 
mates  for  1000  random  samples  of  array  outputs  and  the  accumulation  of  their 
rms  errors.  The  probability  of  detecting  both  emitters  with  the  optimal  de¬ 
tector  was  also  accumulated.  This  was  done  for  .1,  .03,  and  .01  beamwldths 
emitter  separations  and  a  range  of  signal-to-noise  ratios  (SNRs).  The  ac¬ 
curacy  of  the  computations  was  assured  through  a  simple  finite  grid  approxima¬ 
tion  for  the  estimates,  with  no  convergence  problems,  and  through  the  evalua¬ 
tion  of  statistical  confidence  Intervals  for  the -Monte  Carlo  data. 

The  results  of  the  evaluation  indicated  that  the  Cramer-Rao  bound  was 
achievable  by  both  the  MAP  and  WISE  estimators  over  a  wide  range  of  SNR  pro¬ 
vided  a  few  as  10  looks  had  been  taken.  In  general,  the  bound  was  achieved 
wherever  both  signals  were  detectable.  These  results  were  surprising  since 
the  bound  exhibited  unusual  behavior;  for  example,  in  one  SNR  region,  the 
bound  showed  smaller  rms  errors  for  more  closely-spaced  emitters. 

Additional  results  Included  properties  of  the  aposteriori  probability 
density  and  an  analytical  computation  of  the  performance  of  the  known  angles- 
of-arrlval  optimal  detector. 
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The  estimation  of  the  angles-of-arrival  of  wavefronts  incident  upon  a 
sensor  array  is  a  well-known  and  important  problem,  occuring  in  such  varied 
fields  as  geophysics,  oceanography,  emitter  location,  etc.  For  small  angular 
separations  of  the  waves  or  for  small  array  apertures,  estimation  is 
difficult.  As  more  and  more  performance  is  desired  of  such  estimators,  as  the 
cost  of  physically  large  sensor  arrays  increases,  and  as  the  cost  of  real-time 
digital  computation  drops,  the  possibility  of  statistically  optimum  estimation 
becomes  increasingly  interesting  and  practical. 

The  goal  of  this  report  is  not  to  develop  a  practical  optimal  estimator 
but  rather  to  compute  the  performance  of  the  optimal  estimators  and  compare 
this  performance  to  the  Cramer-Rao  bound.  More  specifically,  assuming  a 
Gaussian  signal  model  with  two  Incident  waves  of  known  powers,  the  rras  error 
of  the  Bayes-optimal  MAP  (maximum  aposteriori)  and  MMSG  (minimum  mean  square 
error)  estimators  will  be  compared  via  a  Monte  Carlo  simulation,  and  coiqpared 
to  the  known  powers  Cramer-Rao  bound.  These  rms  errors  will  provide  an 
absolute  lower  bound  on  the  attainable  rms  error  in  angle-of-arrlval  eslmation 
and  will  determine  the  tightness  of  the  Cramer-Rao  bound  for  this  problem. 

1.1  AOA  Estimation  with  Sensor  Arrays 

The  angle-of-arrival  (AOA)  estimation  problem  in  its  most  general  form 
consists  of  the  estimation  of  the  wavevectors  k^  of  say  D  incident  waves 
based  on  the  outputs  of  N  arbitrarily  placed  wavefield  sensors.  Each  sensor 
output  is  some  known  function  of  the  field  in  its  vicinity. 

In  this  report,  we  specialize  the  problem  to  a  simpler  though  common 
case.  Later,  we  assume  only  two  waves  (D“2)  in  order  to  make  computations 
feasible,  but  will  keep  D  arbitrary  whenever  possible.  Assume  that  all  of  the 
T)  waves  are  monochromatic  plane  waves  with  wave-vectors  k^,  complex 
amplitudes  a^  at  the  origin  (r_  *  0),  and  frequency  w0.  The  total  field  at 
position  £  is  then 
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by  superposition  (linear  medium).  Assume  that  the  sensors  are  Isotropic  with 
coherent  detectors  so  that  the  output  of  a  sensor  at  location  jr  Is 


x(r)  -  l  a.e 
1-1 


1*1  *  £ 


(1.2) 


and  that  the  sensors  are  equally  spaced  along  the  x  axis  with  spacing  X/2,  X  - 
wavelength,  so  that  the  position  jr,,  of  the  nth  sensor  Is 


in 


,  n  -  l. 


N 


(1.3) 


Finally,  assume  that  all  the  wavevectors  lie  In  the  upper  x-y  half-plane, 


thus 
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where  Is  the  angle  of  arrival  of  the  1th  wave  with  respect  to  the  y  axis 
(see  Fig.  1.1).  This  Is  necessary  to  make  the  array  output  unique  for  unique 
directions.  Arbitrary  wavevectors  could  be  estimated  by  merely  replacing  the 
array  long  the  y  and  z  axes,  however. 

Plugging  these  assumptions  In  Eq.  (1.2),  we  have  that  the  output  of 
sensor  n  is 
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We  can  rewrite  (1.4)  concisely  by  using  vector  notation.  Define 


"received  signal  vector 


"signal  amplitudes  at  origin  vector 


"direction  vector" 


sin  1-1,  .  .  .  D  "directions-of-arrival" 


I  m  fv(6,)  ...  v(0^)l  "direction  matrix. 


Then  (1.4)  becomes 


x  ■  Va  (1.5) 

We  see  Chat  jt  is  a  linear  combination  of  "discrete-space"  sinusoids,  with 
spatial  frequencies  f^  ■  1/2  0^  *  1/2  sin  Thus,  this  angle 

estimation  problem  is  equivalent  to  spectral  estimation  for  uniformly  sampled 
time  series,  to  within  the  non-linear  transformation  f  *  1/2  sin4>.  In  this 
report,  for  the  sake  of  mathematical  simplicity,  we  will  estimate  the  O^'s 
rather  than  the  <j>i's,  thus  the  relation  to  time  series  frequency  is  a  linear 
one,  f  ■  0/2.  Our  results  are  therefore  directly  applicable  to  time  series 
spectral  estimation  (for  one  look,  see  next  section). 

1.2  Statistical  Model 

Any  real  sensor  array  is  corrupted  by  noise.  Thermal  noise  in  the 
receivers,  for  example,  is  Inevitable.  The  sources,  additionally,  are  best 
modelled  as  noisy  due  to  propagation  environment  effects,  etc.  For  the  usual 
physical  and  analytical  reasons,  these  noises  will  be  taken  as  Gaussian, 
specifically  zero-mean  circular  complex  Gaussian*. 

Accordingly,  our  model  for  the  actual  array  output  is  now 

X  -  Va  +  n  (1.6) 

*z  is  circular  complex  Gaussian,  CN  (m,  o^),  if  Re  z  and  Im  z  are 
jointly  Gaussian  with 

E  (Re  z)  *  Re  m 

E  (Im  z)  "  Im  m 

cov  (Re  z,  Im  z) 


where 


complex  Gaussian  sensor  noise  with 


E(n)  *  0 


zero  mean, 


E(n  nH)  *  I 


unit  power,  and  independent  from  sensor  to 
sensor,  and 


complex  Gaussian  signal  amplitudes  at 
origin  with 


E(a)  -  0 
E(a  aH)  -  P 


zero  mean  and  covariance 

the  ”signal-in-space”  covariance  matrix. 


He  assume  also  that  the  sensor  noise  is  Independent  of  the  signal  amplitude, 
1 .  e  • , 

E(a  nH)  -  [0] 


which  implies  that  the  received  signal  covariance  R  is 
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R  -  E(x  xH) 


E(Va  aH  VH)  +  E(V  a  nH)  +  E(n  aH  VH)  +  E(n  nH) 


V  E(a  aH)  VH  +  0  +  0  +  E(n  nH) 


VPVH  +  I 


(1.7) 


We  also  assume  that  by  sampling  the  array  output  at  a  sufficiently  large 
time  interval  (greater  than  the  slgnal-ln-space's  (finite)  correlation  time), 
we  have  access  to  L  statistically  independent  samples  of  the  random  vector  jc. 


fei J  i-1 


,  L  =  number  of  "looks"  at  array. 


Note  that  this  is  where  the  array  processing  problem  differs  from  the  time 
series  problem,  where  L  is  one  though  N  may  be  large. 

The  estimation  problem  is  to  determine  the  true  directions  of  arrival 
_9*,  given  the  observations  X  and  knowing  apriori  the  powers  P,  the  number  of 
signals  D,  the  number  of  looks  L,  and  the  array  direction  vectors  v(0)  for  all 
9.  This  is  a  classical  parameters-in-the-covariance,  zero  mean  complex 
Gaussian  estimation  problem,  described,  for  example,  in  [15]. 

For  such  a  problem,  there  is  a  simple  sufficient  statistic,  the  sample 
covariance  matrix,  which  obviates  the  need  to  save  all  the  data  vectors  x< , 
i*l,  .  .  .  L.  This  can  be  seen  from  the  probability  densities  (pdf's):  since 
x  is  complex  Gaussian  with  zero  mean  and  covariance  R,  its  pdf  is 
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Hp-i 

-x  R  x 
e  —  — 


Since  the  xj  vectors  are  all  independent  and  identically  distributed,  their 
joint  pdf  is 


p(X)  -  n 


H  D-1 
-x .  R  x. 
— 1  —1 
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where  the  last  step  follow  through  the  use  of  the  trace  Identity 


tr  (AB)  -  tr  (BA) 


Thus,  If  we  define 


.‘4  1  r  H 
R  *  L  l  ii* ! 
1-1 


to  be  the  sample  covariance  matrix,  then 


_  1  -L  tr  R”1  R 

P(X)  -  - - —  e 

irR 


(1.8) 


which  depends  on  the  data  only  through  R.  Thus,  R  is  a  sufficient  statistic 
for  the  collection  X.  This  fact  is  especially  useful  for  Monte  Carlo 
simulations  since  random  samples  of  R  can  be  generated  directly  with  no  need 
to  generate  the  potentially  large  amount  of  data  (L  >  1000  for  some 
applications)  in  X.  All  these  facts  are  well  known  [15]. 

By  way  of  notation,  we  will  often  write  the  expression  p(X)  equivalently 
as  p(x|jB)  or  p(R|jO. 

1.3  Cramer-Rao  Bound 

An  effort  to  determine  the  best  that  one  can  estimate  the  parameters  _0* 
leads  naturally  to  the  Cramer-Rao  bound.  The  Cramer-Rao  bound  [15]  is  a 
well-known  lower  bound  on  the  mean-square  error  of  any  unbiased  estimator. 
Here,  we  assume  that  _6*  is  a  non-random,  though  unknown,  vector  parameter, 
which  leads  to  a  more  useful  bound,  although  the  random  _9*  case  will  be 
considered  in  following  sections. 

The  bound  is  expressed  in  terms  of  the  log  likelihood  function  A(6), 

where 

A(j»  -  In  p(x|j3) 

and  the  Fisher  information  matrix  F,  where 

9A<_0)  3A(j9)  “I 

Fij  "  Ep(x|j})  *  "TeJ"  J  »  1  -  i>  J  -  D 

A 

Then  if  j)  is  any  unbiased  estimate  of  0*,  i.e.. 


the  mean  square  error  of  the  estimate ,  Ag,  given  by 


Ae  ~  E«i  -  I*)  (i-i  )T) 

is  bounded  below  through 
Ae  ^  F'1 


where  the  matrix  inequality  means  that  Ag  -  F"1  is  positive  semidefinite. 
In  particular,  since  a  positive  semidefinite  matrix  must  have  non-negative 
diagonal  elements,  this  implies 


2  A 

o,  - 


> 


For  our  estimation  problem,  A  is  given  from  Eq.  (1.8)  by 

A  »  In  p(x|j»  -  -L  (In |  nR|  +  tr  R-1  R) 

After  much  algebra,  detailed  in  [6],  the  bound  emerges  in  the  form 

A  >  F  *  where 
e  — 


F  -  2  L  Re  {(P-0)  *  (tf*V  -  vfHQW)T  +  (*J  x  (QW)T| 


(1.9) 


A  x  B  ■  matrix  direct  product,  (A  x  B)^  *  A^  B^ 

W  -  vS 


i  i 

V  -  (v  (Oj)  •  ...  i  v(fi^)) 


v(6) 


A 
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3v(e) 

"Te 


.  '  ' 
v  -  (£(e  )  i  ...  i  v(6l)) 

i  i 

w  & 

0  «  (W  +  p"1) 

Notice  the  appearance  of  v( 8) ,  which  is  the  key  dependence  of  the  observed 
data  on  small  changes  in  the  desired  parameter  6,  i.e.,  for  small  emitter 
separations.  Also  notice  that  the  number  of  looks  L  enters  only  as  a  scale 
factor. 

Due  to  the  complexity  of  (1.9),  the  bounds'  behavior  was  examined 
numerically  to  generate  the  curves  in  Fig.  1.2.  Here,  we  assume  the  case  for 
which  performance  evaluations  were  made  in  Chapter  3,  namely  two  uncorrelated 
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Fig.  1.2.  Known  powers  Cramer- Rao  bound. 


signals  of  equal  power.  The  rms  error  bound  on  6},  in  beamwidths  (one 
beamwidth  ■  2/N),  is  plotted  versus  the  array  signal-to-noise  (SNR)  Np,  where 


P  - 


fP 

'•0 


0 

p 


) 


for  various  values  of  the  emitter  separation  A0, 

*  * 

A6  -  0j  -  e2 

This  is  sufficient  since,  by  symmetry,  the  bound  on  0j  equals  the  bound  on 
02,  and  since  by  the  isotropy  of  the  sensors,  the  bound  does  not  depend  on 
the  absolute  position  of  the  emitters  0}  and  02,  but  only  on  their 
separation,  A0. 

We  see  that  the  bound  has  an  unusual  behavior.  In  region  2  of  SNR,  as 
shown  in  Fig.  1.2,  the  rms  error  is  larger  for  larger  emitter  separations, 

which  is  saying  that  widely  spaced  emitters  are  more  difficult  to  estimate 

than  closely  spaced.  This  is  counter  to  intuition  and  also  in  contrast  to  the 
bound  in  the  low  SNR  and  high  SNR  regions  1  and  3,  respectively. 

Also  in  region  2,  we  see  a  plateau  in  the  bound.  This  is  saying  that 
once  a  certain  rms  error  is  attained  (namely,  oj  «  A0/3  at  L  ■  10  looks),  a 
very  large  step  in  SNR  is  required  to  decrease  the  error  further.  Both  of 
these  effects  would  be  of  great  concern  in  the  AOA  estimation  problem  if  they 
were  "real."  But  the  curves  are  only  lower  bounds  on  the  rms  error,  hence  if 
they  are  not  tight  bounds,  the  actual  achievable  rms  error  could  behave 
differently.  For  example,  since  it  is  known  that  for  high  SNR  the  Cramer-Rao 

bound  is  always  achievable  [151,  the  actual  error  might  follow  straight  lines 

coinciding  with  the  bounds'  high  SNR  asymptotes. 

Additionally,  the  curves  are  bounds  only  on  the  etlmation  error  and 
presume  that  it  is  known  aprlori  that  there  are  two  signals  present,  or  that 
two  signals  have  been  detected.  Where  the  rms  error  predicted  by  the  bound  is 


much  larger  than  the  separation  (region  1  in  the  figure),  one  would  not  expect 
both  signals  to  be  detectable. 

These  questions  can  be  resolved  by  computing  the  actual  achievable  rms 
error  and  the  actual  achievable  signal  detectability  (probability  of 
detection),  namely,  the  performances  of  the  optimal  estimator  and  detector. 

1.4  Optimal  Estimation  and  Detection 
1.4.1  Estimation 

The  obvious  way  to  verify  the  tightness  of  the  Cramer-Rao  bounds  of  the 

previous  section  are  through  the  evaluation  of  the  rms  errors  of  the 

statistically  optimal  estimators  ML  (maximum  likelihood),  MAP  (maximum 
apo8teriorl) ,  and  MMSE  (minimum  mean  square  error). 

The  ML  estimate  is  defined  to  be  the  0  which  maximizes  the 

likelihood  function  X: 

A 

%L  *  '  argmax  A(_0) 

_0 

It  has  the  important  property  that  if  any  unbiased  estimate  achieves  the 
Cramer-Rao  bound,  then  the  ML  estimate  does  if  it  is  also  unbiased. 

To  define  the  MAP  and  MMSE  estimators,  we  must  assume  that  6*  is  a 
random  unknown  parameter  with  some  probability  density  p(j)) .  Here,  we  choose 
p(_0)  to  be  a  uniform  density  over  all  possible  _0  values.  Since  0*  - 
sin$i,  we  have  “1  .<  <,  1  i*l,2,  and  since  the  labelling  of  "1"  and  M2” 

is  arbitrary  because  the  signals  are  symmetric  (we  assume  equal  power  here), 
we  also  have,  without  loss  of  generality,  that  02  0j.  Thus,  p(j))  is  a 

uniform  density  on  the  triangle  with  p(j3)  -  1/2  (on  the  left  in  Fig.  1.3). 

Notice  that  this  does  imply  a  nonuniform  density  for  6  »  sin”l0,  which 
is  seen  by  the  change  of  variables  formula  from  probability  theory  to  be 


This  does  not  concern  us  here  but  would  affect 


(on  right  in  Fig.  1.3). 
estimators  of  directly. 

With  p(8)  in  hand,  we  now  have  the  aposteriori  density  of  _0,  p(j)|fc), 
which  is  given  by  Bayes'  rule: 


,  „  p(r|_0)  p(e) 

(6R)  -  — - 


Now  the  MAP  estimator  J^ap  Is  the  maximum  of  this  density 


-MAP  “  arRmax  P(_9|R) 
9 


Because  of  our  choice  of  a  uniform  prior  on  j3,  the  MAP  estimate  is  Identical 
to  the  ML  estimat*. 

Aiap  "  -5ml 

so  we  will,  for  convenience,  consider  only  J^AP* 

The  MMSE  estimate  6mmse  is  the  mean  of  the  aposteriori  density, 

-i»SE  4  •  //  «ep(«|£) 

and  has  the  property  that  It  minimizes  the  total  squared  error,  averaged  over 
the  aposteriori  density  p(_0),  i.e., 

min  E(tr  A  )  ■  min  /  /  dj  tr  A  p(_6) 


iSV'vIvi 


This  Is  not  Che  mean  square  error  Ag  defined  previously,  since  Ag  did  not 
Include  the  average  over  the  prior  p(j>).  The  MMSE  estimator  is  still 
well-defined,  however,  and  could  be  better  than  MAP. 

The  reader  may  be  disturbed  by  our  changing  from  nonrandom  to  random 
J3*.  In  all  of  the  work  which  follows,  _6*  is  actually  assumed  fixed  and 
nonrandom  (though  unknown  to  the  estimator)  because  to  average  the  statistics 
over  many  samples  of  _0*  from  the  density  p(j9)  would  wash  out  effects  of 
Interest,  for  example,  the  dependence  of  the  error  on  _0* .  However,  assuming 
random  j)*  allows  the  definition  and  computation  of  the  conceptually 
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R),  which  describes  exactly  what  is  known  about  JJ 
given  the  observation  R.  By  choosing  a  uniform  prior,  we  make  the  random  and 
nonrandom  _0*  viewpoints  coincide  (e.g.,  *  j^iAp)* 

1.4.2  Detection 

In  order  to  determine  the  regions  where  the  Cramer-Rao  bound  is 
meaningful,  it  is  of  interest  to  examine  the  regions  where  the  two  signals  can 
be  detected.  We  assume  that  the  presence  of  some  signal  is  easy  to  detect  so 
that  the  problem  is  to  distinguish  one  signal  from  two  closely  spaced  ones. 

Let  p(D)  be  the  apriori  probability  of  the  number  of  signals;  thenthe 
minimum  probability  of  error  estimate  of  the  number  of  signals,  D,  is 


argmax  p(D  R) 
D  ! 


where  p(d|r)  is  the  aposteriorl  density  of  D  given  R,  derived  from  Bayes'  rule 


p(RjD)  p(D) 
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and  the  probability  of  error  (Pg)  is 


PE  -  prob  (D  *  D*) 

where  D*  is  the  true  number  of  signals. 

This  detector  is  closely  related  to  the  optimal  estimators,  as  will  be 
seen  in  Chapter  2. 

1.5  Previous  Work  -  Present  Contribution 

Research  relating  to  optimal  angles-of-arrival  estimation  has  been 
carried  out  since  the  early  1960's.  None  of  the  published  results  solve  the 
problem  for  our  particular  combination  of  assumptions,  however.  Most  of  the 
work  assumed  the  easier  deterministic  model,  which  implied  also  that  just  one 
look  was  considered.  Sklar  and  Schweppe  [12]  and  Pollon  [10]  computed  the 
Cramer-Rao  bounds  for  this  problem,  for  various  array  configurations  and  two 
signals.  Young  [17]  assumed  a  prior  density  on  6  as  we  have  and  derived  the 
form  of  an  analog  signal  processor  for  the  MMSE  estimator.  Kslenskl  and 
McGhee  [8]  considered  the  ML  estimator  for  two  signals  and  computed  it 
approximately  by  maximization  over  a  finite  10  by  10  grid  of  directions. 
Pollon  and  Lank  [11]  ran  an  analog  simulation  of  the  two  signal  ML  estimator 
for  a  certain  circular  array. 

Gallop  and  Nolte  [4]  and  Hodgkiss  and  Nolte  [7]  examined  the  one-signal 
case  only,  considered  mainly  the  detection  problem.  Rife  and  Boorstyn  [13, 
14]  simulated  the  one-  and  two-signal  estimators  and  compared  them  to  the 
Cramer-Rao  bound,  but  approximated  ML  by  a  simple  discrete  Fourier  transform. 

The  multiple  look,  Gaussian  model  problem  is  the  more  difficult  one  and 
has  been  little  examined.  El-Behery  and  Macphle  [3]  investigated  the 
one-signal  ML  estimator,  which  reduces  to  a  discrete  Fourier  transform.  They 
performed  an  extensive  Monte  Carlo  Simulation  as  we  will  here  and  found 
attainment  of  the  Cramer-Rao  over  a  surprising  range  of  SNR.  Lalnlotis  [9] 


considered  the  Gaussian  model  for  Bayes'  estimation  in  a  control  theory 

context.  The  estimator  was  computed  by  a  maximization  over  only  10  grid 

points,  however,  one  of  which  was  the  true  parameter  value.  Such  a  choice 

probably  explains  the  performance  observed.  A  control  theory  or  state 
variable  approach  with  finite  grid  approximations  is  also  considered  by  Bucy 
and  Senne  [2],  in  the  context  of  the  phase  demodulation  problem,  and  with  a 
thorough  Monte  Carlo  performance  evaluation. 

The  contribution  of  this  report  consists  principally  of  the  accurate 
computation  of  performance  of  the  MAP  and  MMSE  estimators  for  Gaussian 

two-signal,  known  powers  problem.  The  calculation  is  accurate  enough  to 
distinguish  the  unusual  features  of  the  Cramer-Rao  bound  mentioned  earlier. 

Additionally,  the  performance  of  the  optimal  detector  is  computed  and 
compared  to  an  ad-hoc  approximation  based  on  the  Cramer-Rao  bound,  and  to  a 
rigorous  known-directions  bound. 

In  Chapter  2,  a  more  detailed  derivation  of  the  form  of  the  MAP  and  MMSE 
estimators  and  the  optimal  detetor  for  our  problem  is  given.  The  chapter 

closed  with  some  example  plots  of  p(J)jR)  for  two  emitters.  These  plots 
provide  insight  into  the  nature  of  the  estimation  errors. 

In  Chapter  3,  the  results  of  the  performance  evaluation  are  presented. 
First,  error  sources  in  the  computer  calculation  are  quantified  and  seen  to  be 
acceptable.  Then  the  results  of  the  calculation  are  plotted  and  compared  to 
the  Cramer-Rao  bound  and  the  Pq  bounds.  The  bounds'  behavior  is  seen  to 
indeed  be  real. 

In  Chapter  4,  the  consequences  of  the  tightness  of  the  bound  are 
discussed,  along  with  suggestions  for  further  research. 

The  appendices  contain  derivations  of  the  confidence  intervals  for  the 
Chapter  3  results  and  of  the  known-directions  Pq  bound. 


2.0  THE  MAP  AND  MMSE  ESTIMATORS  AND  THE  OPTIMAL  DETECTOR 

In  this  chapter,  the  equations  for  the  MAP  and  MMSE  estimators  and  the 
optimal  detector  are  derived  in  detail.  Although  the  derivations  are  classi¬ 
cal  [15],  some  simplif ications  specific  to  the  two-emitter  problem  are  ap¬ 
plied.  The  equations  center  on  the  expression  for  the  aposteriori  density, 
P(0.|r),  and  the  chapter  is  concluded  with  example  plots  of  p(j^|R)  illustrat¬ 
ing  some  properties  which  have  been  derived  or  observed. 
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R)  embodies  all  the  information  about  Q_  in  R,  it  is  an  import¬ 
ant  function  itself  in  the  angle  estimation  problem. 

2.1  MAP  and  MMSE  Estimators 

As  mentioned  in  Chapter  1,  the  MAP  and  MMSE  estimates  of  the  direction 
parameters  are  given  in  terms  of  the  aposteriori  density  p(ejft)  by 


4lAP 


argmax  p(e|R) 
6  i 


i«MSE  *  /  /  4  P<«|“) 

The  form  of  p(0^  R)  now  follows  our  Gaussian  assumptions.  We  have  first 
(from  Eq.  (1.8))  that  p(£|§_)  is  a  multivariate  Gaussian  density  with  zero  mean 
and  covariance  R(9): 


p(r|o) 


rR(e)| "L  e"L  tr  1r 


To  define  the  aposteriori  density,  we  oust  use  the  prior  density  of  8^,  p(jO, 
defined  in  Chapter  1.  Then,  by  Bayes'  rule. 


The  notation  here  is  loose  since  this  is  actually  p(x|jO  is  expressed  in 
terras  of  R.  R  itself  is  Complex  Wishart  [5].  '  This  detail  does  not 

effect  p(e|R). 
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Since  p(R)  is  independent  of  6_,  it  is  just  a  normalizing  constant  for  the 
density.  Likewise,  we  chose  a  uniform  prior  density  on  8_  (in  Chapter  1),  so 
that  p(8)  Is  also  a  constant.  Therefore, 


p(£  R)  -  K  p(R|£) 


where  K  depends  on  R  but  not  on  8_. 

We  now  absorb  more  constants  by  using  a  well-known  matrix  inversion 
identity, 


-1 


RCeT1  =  (i  +  VPV11]  -  I  -  vqvh 


where  Q  -  (p  *  +  w)  * 


W  -  vllv 

(same  definitions  as  used  for  Cramer-Rao  bound  in  Section  1.3),  and  get 


:.HV  “ 


p(8ji)  -  K  |»R(.)|-1  .-«•  «  U  -  VQV  )  R 
-  k  |r(9)|-l  rLcr  (R-VQV11*) 


K'  R((0  e 


h: 


-L  L  tr  QV  RV 


where  K1  *  K  i 


-NL  -L  tr  R 


is  once  again  independent  of  8 . 


Finally,  the  determinant  can  be  reduced  since  we  assume  only  two  emitters 


by 


R(8) 


VPVH  +  I 


(1  +  xx)  (1  +  x2) 


1  +  Aj  +  ^  ^  ^2 
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where  and  A2  are  the  two  non-zero  eigenvalues  of  VPVH.  By  a  theorem 
from  linear  algebra,  Aj  and  A2  are  also  the  only  two  eigenvalues  of 
PVHV.  Thus, 


Ai  +  A2  =  tr(PVHV)  -  tr(PW) 


and 

X1  X2  ’  |PVHV|  ’  |™| 


implying 

|R(6_)  j  ■  1  +  tr(PW)  +  |pw| 

which  is  convenient  since  PW  is  only  a  2  by  2  matrix.  Therefore,  we  now  have 


K,  eL  tr  Qv“rV 
(1  +  tr  (PW)  +  |PW|)L 


(2.1) 


This  expression,  unfortunately,  cannot  be  solved  analytically  for  its 
maximum  or  mean,  so  numerical  procedures  must  be  used  to  compute  JmaP  and 
For  our  two-emitter  known  power  case,  0^  is  only  two-dimensional, 
hence,  such  numerical  procedures  are  reasonable.  Of  course,  this  was  the 
principal  reason  for  limiting  the  investigation  to  this  case.  The  procedures 
used  are  described  in  Chapter  3. 


2.2  Optimal  Detector 

The  optimal  detector,  as  also  covered  in  the  first  chapter,  is  given  by 


and  is  the  minimum  probability-of-error  (Pg)  estimate  of  the  number  of 
signals  present.  Since  D  must  be  determined  without  knowing  the  signal 
directions,  this  is  the  composite  hypothesis  testing  problem  as  defined  in 
[15],  with  0.  as  the  nuisance  parameter.  The  test  is  seen  to  be  simply 
related  to  p(6_jR)  of  the  previous  section. 

We  assume  as  in  Chapter  1  that  the  presence  of  some  signal  is  easy  to 
detect,  therefore  we  need  only  distinguish  D-l  and  D*2.  Also,  we  assume  that 
these  cases  are  equally  likely  apriori,  thus  the  aprlori  density  of  D  is 
p(D*l)  *  p(D*2)  ■  1/2.  Finally,  we  must  specify  the  known  power  if  one  signal 
is  present.  The  physically  meaningful  choice  for  this  power  is  pi  +  P2 
where  pj  and  P2  are  the  powers  when  two  signals  are  present.  This  choice 
also  should  be  the  most  difficult  detection  problem  since  otherwise  one  and 
two  signals  could  be  distinguished  on  the  basis  of  total  signal  power  alone. 

Now ,  by  Bayes 1  rule , 


p(R|D)  p(D) 

p(R) 


D  -  1,2 


-  K  P(r|d) 

since  p(D)  and  p(ft)  are  each  independent  of  D. 
unrelated  to  the  K  of  the  previous  section. 

Now,  p(RjD)  is  obtained  from  p(r|d,  B)  by 

p(r|d)  -  /  dQ  p(R|D,  £)  p(0_)  d£ 


The  constant  K.  here  is 


(2.2) 


where  the  integral  is  one  dimensional  for  D“1  and  two  dimensional  for  D«2. 
For  D**2,  p(R|D,  0_)  is  what  we  actually  meant  by  p(R|£)  in  the  previous 
section,  where  two-dimensionality  was  Implicit.  Furthermore,  p(r|d~1,  6)  Is 
obtained  by  evaluating  the  former  along  the  diagonal  0i“92, 


e 

A 


) 
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which  simply  says  that  one  signal  in  direction  6  with  power  pi  +  P2  is 
Indistinguishable  from  two  signals  at  the  same  direction  0  with  powers  pi 
and  p2«  Also,  we  take  the  aprlori  density  of  6  for  D«l,  p(0),  to  be  uniform 
on  -1  <  0  <  l,  hence  p(0)  ■  1/2.  Therefore,  plugging  the  above  and  Eq.  (2.1) 
into  (2.2),  we  have 


p(d|r) 


k*  /  i,  de  P(6)  |r  (  l  )|-1  .*•  “  «*  (  2  )  "RV(  2  )•  D*1 


(2.3) 


K*  //  d£  p(0_)  R(0 )  e 


-L  L  tr  QV(0)  RV(0 ) 


,  D-2 


where  K'  is  chosen  such  that  p(D  -  ljR)  +  p(D  ■  2|R)  ■  1  and  where  the 
determinant  has  been  left  unexpended  for  conciseness.  Thus,  0  is  1  or  2 
according  to  whether  the  1  dimensional  or  2  dimensional  integral  of  p(0_|R)  is 
larger. 

The  Integrals  are,  as  usual,  intractable  and  must  be  computed 
numerically.  Integration  procedures  used  will  be  described  in  Section  3.2. 

2.3  Examples  of  Properties  of  p(£.j&) 

The  aposteriorl  density  reveals  precisely  how  much  information  about  is 
contained  in  R  and  is  the  function  in  terms  of  which  the  optimal  estimators, 
MAP,  MMSE,  and  detector  can  be  expressed.  Therefore,  before  examining  the 
statistics  of  the  estimators  in  the  next  chapter,  it  is  Interesting  to  see 
just  what  it  looks  like. 

First,  we  see  (Fig.  2.1)  that  for  equal  powers  p}  and  P2,  the  pdf  is 

symmetric  about  the  diagonal  0}  -  @2>  This  clearly  must  be  true,  since 

the  numbering  of  the  angles  as  "1"  or  "2”  is  arbitrary.  The  positions  of 

A  A  A 

0)mAP  and  0{qise  are  also  shown  in  the  figure.  Note  that  0mmSE 

particular  oust  be  computed  by  integrating  only  over  the  triangle,  since  the 
mean  of  the  density  over  the  square  (due  to  the  symmetry)  will  always  be  on 
the  diagonal. 

Another  important  feature  we  see  is  that  extraneous  maxima  tend  to  be  in 
the  form  of  ridges  along  the  lines  0i  *  0*  and  @2  M  02»  i.e.  ,  one 


**,  r 


angle  correct  and  the  other  arbitrary.  The  ridges  tend  to  disappear  from  the 
pdf  as  the  number  of  looks  increases  (Fig.  2.2a)  but,  as  can  be  seen  from  the 
log  pdf,  which  is  the  same  as  the  log-likelihood  function  (Fig.  2.2b),  are 
always  present  to  some  extent.  This  property  could  confuse  a  gradient  search 

A 

type  algorithm  for  finding  the  £map* 

If  pi  *  P2,  there  is  an  ordering  on  the  signals,  say  “1”  is  the 

larger  one.  Hence,  the  pdf  is  not  symmetric  and  should  have  a  peak  only  on 

the  correct  triangle  (Fig.  2.3).  Notice  that  the  peak  is  no  longer 

approximately  circular.  In  fact,  it  is  about  10  times  broader  in  the 
direction  as  the  9^  direction,  which  reflects  the  greater  uncertainty  in  the 
smaller  signal’s  location  &2.  Of  course,  if  pi  is  almost  equal  to  p2, 

the  estimator  will  have  difficulty  distinguishing  which  is  actually  larger, 
hence  there  will  tend  to  be  peaks  on  each  triangle  until  many  looks  are  taken 
(Fig.  2.4). 

Finally,  we  note  that  once  the  aposteriori  density  is  essentially 
unimodel  on  the  triangle,  it  can  be  well  approximated  by  a  Gaussian  mound.  In 
fact,  it  is  close  to  a  Gaussian  with  mean  given  by  the  true  signal  locations 
6*  and  variances  given  by  the  Cramer-Rao  bound  (Fig.  2.5).  Notice  that  the 


3.0  MONTE  CARLO  PERFORMANCE  EVALUATION 

In  this  chapter,  we  present  the  results  of  the  Monte  Carlo  evaluation  of 
the  performance  of  the  MAP  and  MMSE  estimators  (their  rms  errors)  and  of  the 
optimal  detector  (Its  probability  of  detection,  Pq).  The  computational  pro¬ 
cedure  and  Its  sources  of  error  are  first  described.  The  errors  arise  from 
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two  sources,  the  piecewise  constant  approximation  to  p(6_|K)  and  the  finite 
number  of  Monte  Carlo  trials  performed.  The  results  of  the  evaluation  are 
presented  next  and  compared  to  the  Cramer-Rao  bound.  The  main  result  is  that 
the  rms  error  of  the  estimators  achieves  the  Cramer-Rao  bound  wherever  the 
Pq  is  usefully  high  and  deviates  from  the  bound  only  for  low  Pq,  where  the 
estimates  become  biased.  This  result  is  surprising  due  to  the  bound's  unusual 
behavior. 

3.1  Computational  Procedure 

Simple  approximations  were  used  in  the  computations  of  performance  to 
make  the  program  fast  and  to  make  the  errors  easy  to  understand  and  quantify. 
The  following  sections  describe  first  the  rms  error  calculation  and  then  the 
Pq  and  bias  calculations. 

3.1.1  RMS  Error  Calculation 

A  .  A 

The  conditional  error  covariance  of  an  estimator  £  ■  £(R)  is,  by  defini¬ 
tion, 

Ae  -  E  ((£  -  3*)  (§.  -  3*)T|if)  ■  J  dft  p(R  0*)  (§  -  0*>  (<[  -  6*)T  (3.1) 

where  0_*  is  the  true  theta  value.  The  estimator’s  rms  error  on  the  it^1 
direction  is  the  square  root  of  the  i^  diagonal  element  of  Ae, 


In  order  Co  numerically  conpute  Ae  from  (3.1),  we  mist  make  basically 
two  different  approximations.  One  is  the  approximate  confutation  of  6  given 
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an  accurate  R  and  the  other  is  the  approximate  integration  with  respect  to  R 
over  its  density.  The  approximations  used  and  their  errors  will  now  be  de¬ 
scribed. 

A  A  A 

Given  an  R,  the  two  estimates  an<*  ®MMSE  are  computed  as 

follows.  First,  p(0jk)/K'  (Eq.  2.1)  is  evaluated  over  a  fixed  finite  grid  of 
0  values,  say  G  by  G,  uniformly  spaced  over  the  square  (see  Fig.  3.1).  Since 

I”  ~ 

K)  is  zero  everywhere  except  in  a 

small  region  around  6_*,  the  square  is  taken  to  be  the  smallest  one  which 
safely  encloses  the  pdf's  support.  This  is  actually  (-1,  1J  x  1-1,  1]  only 
for  the  low  SNR  cases  examined.  Then  this  is  numerically  Integrated  over  the 
upper  triangle  to  find  K' .  6^mse  is  now  the  numerical  integral  of  £>p(0_|R) 
over  the  triangle  and  is  the  location  of  the  grid  point  with  the 

largest  value  of  p(£|k).  The  integrals  are  taken  such  that  they  are  exact  if 

I'N  ' 

R)  is  assumed  to  be  piecewise  constant,  i.e.  ,  constant  on  each  grid  cell, 
as  shown  in  Fig.  3.1.  This  implies  that  the  computation  is  a  simple  summation 
plus  a  few  correction  terms  and  that  the  only  error  source  is  the  degree  to 
which  p(ejR)  is  not  piecewise  constant.  The  error  will  go  to  monotonically 
zero  as  G  is  Increased  since  this  approximation  converges  uniformly  on  the 
square  to  p(_0_)  R)  [2]. 

In  fact,  a  plot  was  made  of  the  integrals  versus  G  to  determine  how  fine 
a  grid  was  required  (Fig.  3.2).  For  typically  smooth  true  pdfs,  G  50  was 

A 

adequate  for  an  accuracy  of  about  2%  of  the  true  value.  Errors  in  J^oise  an<* 
were  not  surprisingly  of  the  order  of  the  grid  cell  size.  Making  these 
errors  about  2%  places  them  below  the  errors  due  to  the  second  source,  the  in- 

A 

tegration  with  respect  to  R. 

This  integration,  involving  as  it  does  the  complicated  function 
is  hopeless  to  perform  analytically.  It  is  also  of  such  dimension¬ 
ality  (e.g.,  81-dimensional  for  nine  antennas)  that  the  only  tractable 
numerical  integration  method  is  the  Monte  Carlo  method  [16]. 


GRID  POINT 


REGION  OF 
INTEGRATION 


APPROXIMATE  pdf  =  TRUE  pdf  AT  GRID  POINTS 


Fig.  3.1.  Piecewise-constant  approximation. 


Monte  Carlo  integration  in  this  case  amounts  to  the  simulation  of  the 
estimators,  computation  of  the  errors  and  squared  errors,  and  accumulation  of 
their  sample  means,  i.e. , 

A  -  -  l  (§(R.)  -  0*)  (6 (R  )  -  e*)T 

e  M  W  -  1  “  “  1 

A 

where  M  is  the  number  of  statistically  independent  Monte  Carlo  trials,  R^  is 
the  sample  covariance  chosen  from  the  population  p(R|£_*)  on  the  ifch 
trial,  and  Ae  is  the  estimated  error  covariance. 

/s 

Now  Ae  is  a  random  variable  with  some  unknown  probability  density.  Its 
expected  value  is  the  true  Ae  and  its  variance  is  determined  by  its  pdf.  If 
we  assume  that  its  pdf  is  Gaussian,  which  will  be  approximately  true  for  large 
M,  then  its  variance  can  be  computed  and  used  to  compute  the  number  of  trials 
required  for  a  specified  accuracy.  Choosing  a  99%  confidence  interval  width 
of  ±2  dB  as  adequate  for  distinguishing  the  interesting  features  of  the 
bounds,  we  found  by  making  trial  runs  that  about  1000  Monte  Carlo  trials  were 
required  (see  Appendix  A). 

Since  the  rms  errors  for  each  direction  6^  and  0  2  are  equal  by  symmetry, 
the  final  estimate  of  the  rms  error  for  either  direction  was  taken  to  be 


A  A 


where 


3.1.2  Probability  of  Detection  Computation 
The  probability  of  detection  (Pp)  for  the  optimal  detector  is  given  by 


y.< 


where 


£  -  /*  Z1  dflj  ae2  p(e_|R)  -  /*  de  p((J)|r) 


(3.3) 


and  p(£  |£  )  is  £'s  probability  density.  So  Pq  is  the  probability  that  the 
integral  over  the  triangle  is  greater  than  the  integral  along  the  diagonal. 

Unfortunately,  p(l  |  ®.*)  is  not  tractable  to  compute  analytically  and 
probably  not  easy  to  integrate  as  required  in  (3.2).  So  Pq  must  be  computed 
by  a  Monte  Carlo  simulation*. 

We  define  the  0-1  random  variable 


! 


0,  £  <  0 

1,  £  >  0 


(3.4) 


i.e.,  £  ' 


1  iff  two  signals  are  detected.  Then  Pq  is  the  expected  value  of 


E(£ ' ) 


Therefore,  Pq  can  be  approximated  by 


-  i  M 

Pd'H  w 


£ ' 
i 


(3.5) 


where  M  is  the  number  of  statistically  independent  Monte  Carlo  trials,  as 
before,  and  £^  is  a  sample  of  £'  from  the  population  of  p(£ '  | £*)  on  the 
ith  trial.  £^  is  generated,  of  course,  by  computing  pC®_j R-i )  >  finding  its 
2-D  and  1-D  integrals  as  in  (3.3),  and  using  the  definition  (3.4). 


*lf  it  is  assumed  that  the  detector  knows  6*,  then  this  Pq  can  be  computed 
and  serves  as  a  convenient  upper  bound  for  the  present  one  (see  Appendix  B). 


The  approximations  made  here  are  In  the  piecewise  constant  approximation 
for  p(£  R)  in  the  integrals  and  in  the  finite  variance  of  the  random  variable 
Pq.  Assuming  the  grid  has  already  been  chosen  fine  enough  to  make  the  inte¬ 
grals  for  the  rms  error  calculation  in  the  previous  section  sufficiently  accu¬ 
rate,  only  the  Monte  Carlo  error  remains. 

As  explained  in  Appendix  A,  this  error  is  quantified  by  the  observation 

A 

that  Pq  has  approximately  a  Gaussian  pdf.  Using  the  1000  trials  required 
for  the  rms  error  calculation,  we  find  in  the  Appendix  that  the  99 X  confidence 
interval  for  Pp  is  about  ± .05. 


3.1.3  Bias  Computation 
The  bias  b  of  an  estimate  9  is  defined  as 


b  =  E(0_  -  0) 

and  was  computed  for  MMSE  and  MAP  exactly  the  same  way  as  the  rms  error  and 
the  probability  of  detection.  Given  M  sample  angle  estimates,  the  estimate 
for  the  bias  was 


A  1  V  c A  *•> 

b=M  l  • 


Confidence  Intervals  for  b  follow  easily  and  are  derived  in  Appendix  A. 


3.2  Results 

A  computer  simulation  program  was  written  which  performed  the  Monte  Carlo 
evaluation  of  rms  errors,  probability  of  detection,  and  bias  of  the  optimal 
estimators  (MAP  and  MMSE).  The  results  of  running  that  program  for  various 
array  signal-to-noise  ratios,  emitter  separations,  and  numbers  of  looks  are 
presented  in  this  section.  In  every  case,  two  equal  power,  uncorrelated  emit¬ 
ters  and  a  uniform  linear  aray  with  nine  or  five  elements  were  assumed.  The 
results  indicate  that  rms  errors  and  detection  thresholds  corresponding  to  the 
Cramer-Rao  bound  are  in  fact  achievable  given  as  few  as  10  looks  at  the  array. 


Incidentally,  the  running  time  for  the  program  for  9  elements  and  a  50  by 
50  grid  was  about  0.19  sec/trial  on  an  Amdahl  470  and  3.1  sec/trial  on  a  VAX 
11/750. 

3.2.1  RMS  Error  Results 

Figures  3.3a  through  3.3c  show  the  computed  rms  errors  of  the  MAP  and 
MMSE  estimators  versus  SNR  for  emitter  separations  of  .1,  .01/  lO,  and  .01 
beamwidths,  respectively,  and  all  at  10  looks.  Figure  3.3d  is  an  overlay  of 
3.3a,  b,  and  c.  The  array  for  these  cases  has  nine  elements  implying  a  beam- 
width  of  2/9  m  .222.  For  an  emitter  separation  of  A6,  the  true  signal  loca¬ 
tions  were  taken  to  be 


-TA9 


i.e. ,  centered  on  the  array  boresight. 

Below  about  5  dB  in  SNR,  we  see  from  the  figure  that  the  rms  error  is  in¬ 
dependent  of  the  true  signal  separation.  In  fact,  it  is  independent  of  almost 
everything  in  this  region  since  the  aposteriori  pdf  is  approximately  uniform 
over  the  triangle,  on  the  average.  Implying  that  the  estimates  vary  approxi¬ 
mately  uniformly  over  the  triangle.  An  exactly  uniform  pdf  would  imply  an  rms 
error  of  4.5//X*  2.6  beamwidths,  approximately  that  of  MAP.  MMSE  is  somewhat 
better  since  it  tends  to  choose  estimates  in  the  center  of  the  triangle  given 
a  uniform  pdf ,  and  this  happens  to  be  close  to  where  the  signals  actually 


Above  about  20  dB  is  the  more  interesting  region  of  SNR.  First,  we  see 
that  the  MAP  and  MMSE  estimator  errors  are  identical.  This  is  not  surprising, 
since  the  aposteriori  density  is  substantially  unimodel  and  symmetric  about 
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its  mode  in  this  region,  so  that  its  mean  and  peak  coincide.  More  signifi¬ 
cantly,  we  see  that  the  rms  error  is  always  less  than  or  equal  to  the  Cramer- 
Rao  bound  (when  less,  the  estimator  was  biased)  and  tracks  the  bound  closely 
beyond  its  first  "knee,"  which  we  will  later  see  is  for  (L“10)  at  the  detec¬ 
tion  threshold.  At  40  dB  SNR,  for  example,  we  see  that  the  rms  error  is 
clearly  (including  the  confidence  interval)  smaller  for  the  smaller  separa¬ 
tions,  verifying  this  unusual  behavior  of  the  bound  noted  in  Chapter  1. 
Additionally,  we  see  that  the  rms  error  is  approximately  constant  as  the  SNR 
increases  from,  for  example,  50  to  70  dB  for  the  .01  beamwidths  emitter  sepa¬ 
ration,  verifying  the  plateau  in  the  bound.  Therefore,  these  counterintuitive 
aspects  of  the  known  powers  Cramer-Rao  bound  observed  in  Chapter  1  apply 
equally  well  to  the  behavior  of  the  optimal  estimators.  It  was  not  “just  a 
bound"  after  all. 

Of  course,  the  preceding  observations  are  for  10  looks,  but  this  is  a 
surprisingly  small  number:  the  convergence  of  the  rms  error  to  the  bound  as 
the  number  of  looks  increases  is  quite  rapid.  Figures  3.4a  and  3.4b  show  this 
convergence  for  two  selected  SNRs  and  separations,  50  dB,  .01  beamwidths,  and 
30  dB,  .1  beamwidths,  respectively.  The  bound  is  achieved  to  within  the  99% 
confidence  interval  after  about  10  looks  and  continues  to  track  the  bound  for 
higher  numbers  of  looks. 

Expecting  that  the  number  of  looks  required  to  achieve  the  bound  may  be 
related  to  the  number  of  antennas,  a  run  was  made  at  50  dB,  .01  beamwidths  for 
only  five  antenna  elements  (Fig.  3.5).  As  can  be  seen  from  the  figure,  the 
threshold  is  about  the  same  as  for  nine  antennas.  It  could  be  that  the 
threshold  depends  on  the  number  of  signal  directions  one  is  trying  to  esti¬ 
mate. 

3.2.2  Probability  of  Detection  Results 

Figures  3.6  and  3.7  show  the  probability  of  detection  confuted  versus  SNR 
for  the  same  three  separations  and  the  same  number  of  looks  as  used  previous¬ 
ly.  From  these  plots,  we  see  the  region  of  SNR  for  which  the  rms  errors  are 
of  any  practical  interest. 
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Fig.  3.7.  Comparison  of  optimal  Bayes  detector  performance  to  known-0  bound 


Below  10  dB  In  SNR,  we  see  that  the  Pq  does  not  depend  on  the  true  sep¬ 
aration,  as  was  the  case  for  rms  error.  The  detection  here  is  probably  no 
better  than  a  random  choice. 

At  higher  SNR,  we  see  sharply  defined  detection  thresholds  which,  in 
fact,  occur  at  the  point  where  the  rms  error  began  to  track  the  Cramer-Rao 
bound.  If  we  take  the  threshold  to  be  when  Pq  .5,  then  a  simple  formula 
based  on  the  figures  for  the  threshold  SNR  versus  separation  is 

Detection  Threshold  SNR  «  —  _•  (3  5) 

(A0  ) 

where  Ad  *  emitter  separation  in  beamwldths,  valid  for  L  »  10  looks.  This 
formula  yields  20  dB,  30  dB,  and  40  dB  SNR  for  .1,  .01  / 10,  and  .01  beamwidths 
separation,  respectively. 

The  dotted  curves  in  Fig.  3.6  are  an  ad-hoc  approximation  to  Pq  based 
on  the  Cramer-Rao  bound.  They  were  generated  by  assuming  that  p(6|ft)  was  a 

JL  * 

Gaussian  mound  with  mean  given  by  and  covariance  given  by  the  Cramer-Rao 
bound  for  the  separation  and  SNR  used.  The  integrals  of  the  density  over  the 
upper  triangle  (call  this  12)  and  along  the  diagonal  (II)  were  evaluated 
exactly  as  in  Eq.  (2.3).  The  approximation  to  Pq  was  then  taken  as 


Although  this  derivation  is  loose,  here  we  see  that  it  works  quite  well. 
This  implies  that  the  aposterlori  density  is  fairly  well  approximated  by  a 
Gaussian  mound  with  the  Cramer-Rao  bound  as  its  covariance.  This  is  also  say¬ 
ing  that  the  off-diagonal  term  of  the  Cramer-Rao  bound  is  being  achieved, 
since  it  corresponds  to  the  orientation  of  the  Gaussian  mound,  hence  to  the 
relative  values  of  12  and  II  and  thus  to  the  Pq.  Such  was  ejected  from 
examining  the  picture  of  the  pdf  for  various  cases,  as  in  Fig.  2.5. 


In  Fig.  3.7,  the  confuted  Pq’s  are  compared  to  the  known  0_  Pq  men¬ 
tioned  in  Chapter  2  and  derived  in  Appendix  B.  It  is  indeed  an  upper  bound, 
although  not  a  very  tight  one.  The  known  9_  curves  also  follow  the  1  /  ( A 0 ) ^ 
dependence  in  Eq.  (3.1)  closely,  indicating  that  knowing  0_  is  equivalent  to  a 
fixed  gain  in  SNR,  independent  of  separation. 

In  Fig.  3.8,  Pq  is  plotted  versus  the  number  of  looks,  for  the  same  two 
cases  as  in  Fig.  3.4.  The  threshold  here  is  at  about  five  looks,  which  agrees 
fairly  closely  with  the  threshold  observed  for  the  rms  error  to  achieve  the 
bound. 

3.2.3  Bias  Results 

In  Figs.  3.9a  through  3.9c,  the  bias  of  the  MAP  and  MMSE  estimators  ver¬ 
sus  SNR  is  plotted  as  a  percentage  of  the  rms  error.  Below  »20%,  one  can  con¬ 
sider  the  estimator  unbiased  since  the  bias  affects  the  rms  error  only  by 
(.2)2  »  4%.  From  the  plot,  we  see  that  the  estimator  is  indeed  biased,  as 
it  must  be,  in  the  low  SNR  region  where  its  rms  error  was  below  the  Cramer  Rao 
bound,  but  becomes  unbiased  at  about  the  detection  threshold. 

The  confidence  interval  shown  is  so  wide  because  it  was  taken  as  the 
maximum  (most  conservative)  estimate,  i.e.. 
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Fig.  3.8.  Performance  of  optimal  Bayes  detector  versus  number  of  looks 


4.0  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

The  most  important  result  of  this  research  is  the  presentation  of  the 
agreement  between  the  rms  errors  of  the  optimal  estimators  and  the  known 
powers,  Gaussian  signal  Cramer-Rao  bound.  The  two  unusual  features  of  the 
bound,  that  the  rms  error  is  constant  over  a  range  of  SNR  and  is  smaller  for 
larger  emitter  separations  over  about  the  same  SNR  range,  have  been  verified. 
The  verification  was  made  rigorous  through  the  use  of  a  straightforward 
finite-grid  computation  for  the  MAP  and  MMSE  estimates  and  through  the  evalua¬ 
tion  of  statistical  confidence  intervals  for  Monte  Carlo  results. 

As  a  consequence  of  this  agreement,  the  Cramer-Rao  bound  can  be  now  used 
with  some  confidence  to  predict  the  ultimate  performance  of  angle-of-arrival 
estimators  for  various  combinations  of  assumptions.  This  is  valuable  since 
the  bound  is  considerably  easier  to  confute  than  the  optimal  estimator  per¬ 
formance.  One  application  of  this  was  the  ad-hoc  formula  for  the  Pq  derived 
from  the  Cramer-Rao  bound  in  Chapter  3,  which  agreed  well  with  the  Py  com¬ 
puted  by  the  simulation.  Without  running  additional  simulations,  the  bound 
can  be  used  for  cases  of  different  numbers  of  sensors,  looks,  or  different 
array  configurations  from  those  run  in  this  research  to  predict  performance. 
The  unknown  powers  Cramer-Rao  bound  may  even  be  used,  with  some  confidence,  a 
case  for  which  the  simulation  is  considerably  more  difficult,  requiring  a 
4-dimensional  grid. 

Although  it  was  the  non-random  parameters  Cramer-Rao  bound  which  was 
actually  verified,  a  random  parameter  (Bayesian)  approach  to  the  overall  esti¬ 
mation  problem  was  used  in  order  to  define  the  MMSE  estimator,  for  one,  and 
more  Importantly  the  aposteriori  density  of  the  directions-of-arrival 
p(£|R).  The  properties  of  this  function  are  important  since  they  could  lead 
to  efficient  procedures  for  computing  the  optimal  estimators  and  to  insights 
into  the  sources  of  estimation  errors. 

While  the  known  powers  problem  is  certainly  of  theoretical  interest, 
Isolating  as  it  does  the  difficulty  in  angle  estimation  alone,  it  might  not  be 
of  great  practical  utility.  Further  research  should  therefore  center  on  the 
unknown-powers  problem.  This  research  could  be  approached  in  several  ways. 


A  straightforward  extension  to  a  4-dimensional  grid  search  is,  of  course,  not 
computationally  feasible.  One  could  test  the  importance  of  knowing  the 
powers  by  providing  the  known  powers  estimator  with  incorrect  powers,  which 
would  determine  how  finely  tuned  the  estimator  is  to  the  particular  signal 
model  assumptions.  If  it  is  not  very  sensitive  to  misinformation  about 
powers,  then  one  would  expect  completely  unknown  powers  to  cause  little 
additional  degradation. 

Extensions  of  this  research  to  even  more  general  cases,  such  as  >  2 
signals  or  to  2-dimensional  angles  parameters,  await  sophisticated  but 
provably  accurate  procedures  to  search  a  many-dimensional  space.  Our 
restriction  to  just  two  known  power  signals  was  made  specifically  to  avoid  all 
questions  of  convergence  which  typically  arise  with  such  procedures. 
Hopefully,  the  results  presented  here  will  be  of  some  help  in  the  general 
problem. 

Ia 

R),  for  example,  provide  some  insight  into  the  animal 
which  must  be  searched.  The  fact  that  p(9.|h)  seems  to  approach  a  Gaussian 
mound  with  the  Cramer-Rao  bound  as  its  covariance,  if  quantified,  could  be 
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APPENDIX  A 

CONFIDENCE  INTERVALS  FOR  THE  MONTE  CARLO  RESULTS 


In  this  appendix,  we  derive  the  confidence  Intervals  for  the  estimates  of 
bias,  mean  square  error,  and  probability  of  detection  shown  In  the  figures  of 
Chapter  3.  Much  of  this  appendix  Is  abstracted  from  Chapter  4  of  [1]. 

Let  a  be  the  true  value  of  some  parameter  to  be  estimated  and  let  a  be 
Its  estimate.  Then  a  100y  %  confidence  Interval  for  a  is  a  pair  of  random 
variables  i(o)  and  u(a)  such  that 

prob  (i(aK  a  <  u(o))  ■  y 

[18,  p.  363].  Here  we  choose  the  symmetric  interval  where  £  and  u  are  of  the 
form 

.A  A 

£(a)  *  a  -  k 
u(a)  a  a  +  k  , 

thus , 

prob(  |a  -  aj  <  k)  «  y 

To  compute  k.  It  suffices  to  know  the  pdf  of  z  ■  a  -  a,  say  p(z).  This 
Implies  that  k  Is  given  by 

k 

/  dz  p(z)  ■  y  • 

-k 

In  each  of  our  cases,  we  are  estimating  the  mean  x  of  a  random  variable  x 
given  M  independent  samples  {xi,  .  .  .,  x^}  of  that  variable.  Hence,  the 
best  estimate  of  x  is 
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If  we  assume  now  that  M  is  so  large  that  x  is  approximately  Gaussian, 
x  is  Gaussian  with  mean,  variance  given  by 


then  x  - 


E<x  "  x)  - 
2 

o  “  var(x)  = 


l  E(xt)  -  x 


2  -2 
E(x  J  -  yT 

M 


0 

var(x) 
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and  we  can  easily  derive  a  symmetric  confidence  interval.  In  particular,  the 
99.74%  confidence  interval  is  given  by 


x  -  x  3o 


Unfortunately,  although  the  pdf  of  x  may  be  very  nearly  Gaussian,  the  pdf  of  x 
is  unknown  so  is  also  unknown.  Here,  we  get  an  approximate  confidence 
interval  by  using  the  estimate  of  a ^ 


2 

a 


-2 


x 


which  should  be  a  good  approximation  for  large  N. 

For  bias  b  and  mean  square  error  s,  we  are  estimating  the  means  of 

* 

x  *  e  -  e 

and 

x  «  (0  -  6*)2 

A  A 

respectively.  Hence,  given  M  sample  angle  estimates  {6i,...,6m}>  we 
compute  the  estimates 


6  -  ±  KSj  -  8*) 

5  -S  I(Si  -  9*>2 
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and  the  fourth  moment 


y* 


K 


*  1  r  ✓  £  *  v  4 

u4  "  M  ^(0i  "  6  ) 

implying  that  the  approximate  99.74 X  confidence  intervals  as  derived  above  are 


|f>  -  b|  3^8  -  £2 


The  interval  for  the  rms  error,  /s"!  was  derived  from  the  interval  for  s  as 

Vs  ~  f  V*"”  <VS  +  ^ 

On  the  plots,  all  of  the  intervals  (for  each  SNR  and  separation)  were 
approximately  the  same  so  the  largest  one  was  taken  as  the  interval  for  all 
the  points.  By  making  runs  for  various  values  of  M  and  computing  the 
confidence  intervals,  M-1000  trials  was  found  adequate  for  a  ±2  dB  interval 
for  rms  error. 

For  the  probability  of  detection,  we  can  be  a  little  more  precise  since 
we  are  estimating  the  mean  of  a  binomial  random  variable 

!0,  no  detection,  probability  1  -  Pq 
1,  detection,  probability  Pj) 

and  can  compute  the  worst  case  var(x): 

var(x)  -  E(x2)  -  E(x)2  -  PD  -  P2  -  PpU  -  Pp) 

-  7  *  I  “  k  for  all  0  <  PD  <  1 

implying 

if 


Therefore,  the  worst  case  99.74%  confidence  interval  for  Pq  is  (stil 
assuming  Pq  is  approximately  Gaussian) 


2/M 

which,  at  M  *  1000  is 

|fD  '  PdI  <  -047  • 


as  shown  in  the  Pq  plots  of  Chapter  3. 


APPENDIX  B 

THE  KNOWN  ANGLE -OF -ARRIVAL  DETECTOR 


The  probability  of  detection  (Pp)  of  the  optimal  Bayes  detector 
described  in  Chapter  2  can  be  upper  bounded  by  the  Pq  of  a  detector  which 
knows  a  priori  the  angles  of  arrival  _0.  This  latter  Pp,  unlike  the  former, 
may  be  computed  analytically  and  hence  is  a  convenient  bound.  In  this 
appendix,  we  derive  in  detail  the  P^  of  the  known  _0  bound,  which  is  seen  to 
be  the  solution  of  the  general  Gaussian  hypothesis  testing  problem  in  [15]. 

By  the  phrase  "known  _9,"  we  mean  more  precisely  that  if  there  are 
actually  two  signals  present  (hypothesis  H2),  we  know  their  angles 


and  if  there  is  only  on_  signal  present  (hypothesis  Hj),  we  know  its  angle 
9.  Of  course,  we  assume  throughout  this  report  that  the  powers  of  all  the 
signals  are  also  known,  say  pi  and  p2  for  two  signals  present,  and  p  for 
one  signal.  Normally,  for  the  problem  to  make  physical  sense,  we  assume 


p  “  P1  +  p2 

(B.l) 

Pl°l  +  P202 

9  *  - — -  (center  of  power). 

P1  +  p2 


Knowing  the  angles  and  the  powers  now  means  that  we  completely  know  the 
covariance  of  the  observations  (x^  on  each  hypothesis.  Therefore,  the 

detection  problem  is  to  discern  based  on  the  L  independent  samples  {xjJibI 
whether 


HjS  ~  CN  (0,  Rj) 
or 

H2:  ™  (0,  R2) 

where 

R1  -  1  +  pv(0)  v(0)H 
R2  -  I  +  V(_0)  PV(_9)H 


(B.2) 


_v(  0)  *  direction  vector 

i 

V(0)  -  (v(  0^ )  '  v(02)) 

CN(m,  R)  *  complex  normal  pdf  with  mean  m  and  covariance  R. 

This  Is  the  well-known  equal  mean,  differing  covariance  Gaussian  hypothesis 
testing  problem  [15,  p.  113,  case  1C].  The  computation  of  its  Pj)  is 
tedious,  but  straightforward. 

The  likelihood  ratio  test  is 


/■a  ON 


where  y'  ■  1  for  minimum  probability  error  (Pg>  detection  assuming  each 
hypothesis  is  equally  likely  (P(Hj)  -  p(H2>  “  1/2). 

Plugging  the  complex  normal  distribution  into  (B.3)  and  using  the 
independence  of  the  Xi ' s ,  we  get 


*  Rl|  I  V  <  H  „-l  H  „-l  A  ^ 

— jr  exp  l '  Jml  r2  ^i  •  Ri  it)  1  <  y' 


4  (RI1  ~  *2l)  ^i  <*  to  y'  -  L  to  —■ 


L  tr  ((Rj1  -  R^1)  R)  \  <n  y’  +  L 


or,  finally. 


((Rj1  -  R'1)  R)  \  -55L.il  + 


(B.4) 


where 


R-r  l  xt4 

b  i-i  11 


-  R"1)  R) 

be  the  statistic  and 


+  to 


be  the  threshold.  Then  the  Pj)  of  this  test  is 


PD  =>  prob(i>y)  ■  /  P^CuIh  )  du 

Y 

where  p£  ( •)  is  the  pdf  of  l,  which  will  be  obtained  through  the 
characteristic  function  of  l. 

The  characteristic  function  of  £, 

<t>z  (jw)  - 

is  given  on  hypothesis  by 


♦t  (Hv 


1  -r  "I1 


-L 


as  derived  in  [5,  Lemma  4.1J. 


Pulling  out  an 


*i  <Hv 


i 


L 


Now  let 


UJ\U  -  r£/2  (r"1  -  r"1)  r£/2 


(B.5) 


be  the  eigenvalue  decomposition  of  the  Hermetian  matrix  on  the  right-hand 
side.  Then, 

-L 


h  (Hv "  1 1  ■  uauH| 


I  -  1“  A 
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-L 


1  -i“  X 

1  L  A1 


-L 


det 
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This  is  a  difficult  characteristic  function  to  invert  until  we  note  that  only 
three  eigenvalues  are  non-zero  since 
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rank  (r^  (r”1  -  Rp  R^) 


Now 


*  rank  (R^  -  R^1 ) 

■  rank  (i  -  q  _v  vH  -  (I  -  V  Q  VH)) 

■  rank  (v  0  V**  -  q  v  v" ) 

•  3  at  most. 


«*|V 


(-  1) 


$$$(•-  xr)  (s  -  x2A) 


,-i 


which  may  be  expanded  into  partial  functions 


VL*|V 


1 


3  L 
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where  bj  ■  1/Xj  and  the  coefficients  a^j  are  given  by 
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or,  finally. 
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where 


b3  “  bl  _  (2L  -  i  -  1)!  r2L  -  i  -  1 

r  "  b2  -  ^  ’  uo  "  (L-i)!  (h  ~  1)!  1  L  -  i  J 

(L  +  £  -  1)!  (2L  -  i  -  £  -  1)!  (L  +  £  -  1)  (L-i-£+l) 
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and  ai2  and  ajj  are  obtained  symmetrically. 

Given  these  a^j,  we  now  inverse  transform  Eq.  (B.6): 


P £(*>  -  Ff*£(J 


L 
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Therefore, 
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and  P  ( •,  •)  is  the  Incomplete  Gamma  function 


(B.7) 


(sgn  y) 


(n,  x)  -  /  dt  tn-1  e~l 


Equation  (B.7),  of  course,  gives  the  probability  of  false  alarm  (Pp) 
also  if  the  eigenvalues  Aj  are  those  corresponding  to  -Rj  instead  of  R£ 
in  Eq.  (B.5). 

A  RATFOR  program  was  written  which  computes  Pj)  and  Pp  for  arbitrary 
choices  of  0j,  02,  Pi,  P2t  and  L  from  Eq.  (B.7),  and  was  used  to 

generate  the  solid  curves  in  Fig.  B.l. 
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Fig.  B.l.  Known- 0  probability  of  detection  comparison  of  exact  and 
Gaussian  computations. 


Unfortunately,  there  are  numerical  problems  in  computing  the  a^j's. 
Fortunately,  the  pdf  of  1  approaches  Gaussian  for  large  numbers  of  looks  and 
Is  reasonably  close  to  Gaussian  at  the  point  where  the  exact  computation 
falls.  Accordingly,  the  program  optimally  computes  Pq  and  Pp  based  on  a 
Gaussian  approximation  given  by 
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This  approximation  appears  as  dotted  curves  in  Fig.  B.l 
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